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Abstract 


Through the method of truncation functions, the oceanic geo id undulation 
is divided into two constituents: an inner zone contribution expressed as an in- 
tegral of surface gravity disturbances over a spherical cap, and an outer zone 
contribution derived from a finite set of potential harmonic coefficients. Global, 
average error estimates are formulated for undulation differences, thereby pro- 
viding accuracies for a relative geoid. The error analysis focuses on the outer 
zone contribution for which the potential coefficient errors are modeled on the 
assumption that the coefficients are determined from a global distribution of 
l°x 1° mean anomalies. The method of computing undulations based on gravity 
disturbance data for the inner zone is compared to the similar, conventional 
method which presupposes gravity anomaly data within this zone. The two meth- 
ods exhibit analogous error characteristics, the estimated errors of the gravity 
disturbance method being only slightly better. For continuous and errorless 
gravity data inside a spherical cap having a radius of 10° and with potential 
coefficients derived to degree 180 from a global set of l°x 1° mean anomalies 
given to an accuracy of ±10 mgal (± 1 mgal), the typical error in the difference 
of undulations is 30 cm to 40 cm (5 cm to 10 cm), depending on their separation. 

In the absence of more detailed cap data, the error is as high as 160 cm (70 cm). 

In the latter case, the corresponding error for l°x 1° mean undulation differences 
is about 140 cm (30 cm). 
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1. Introduction 


This report again attacks the problem of determining the geoid accurately 
from gravity information over the entire globe. The route that is normally 
followed starts with a formulation of the geoid undulation in terms of mean 
gravity anomalies, or potential harmonic coefficients, or both and follows with 
an examination of how the errors in these observed quantities propagate into the 
computed undulation; see for example, Rapp and Rummel (1975) and Christo- 
doulidis (1976), to mention only two. 

Christodoulidis based his thorough error analysis on the method by which 
undulations are computed from mean gravity anomalies and a set of potential 
coefficients (according to Molodenskii's truncation theory). Moritz (1974) out- 
lined a procedure in which the undulations are determined not by integrating 
gravity anomalies (Stokes' formula), but by integrating gravity disturbances . 

It is the purpose of this report to carry out an error analysis of undulations 
(more precisely, undulation differences) that are also obtained according to 
Molodenskii's truncation theory, but as applied to an integral of gravity distur- 
bances. The structure of the analysis is very similar to that of Christodoulidis. 
From his conclusions, it is evident that the procurement of a highly accurate 
(~ 10 cm) relative geoid on the oceans requires a dense and accurate network of 
gravity anomalies around each point of computation, as well as potential coeffi- 
cients to at least degree 70. We shall see that the method of determining the 
geoid from gravity disturbance data under most practical circumstances can 
claim little, if any, substantial improvement in accuracy over the conventional 
method. Because the ultimate determination of the sea surface topography is 
foreseen to be effected by combining accurate geoid undulations and satellite 
altimetry, the areas of primary interest are the oceans. Therefore, this al- 
ternate method does offer two advantages; first, measured gravity values on 
the ocean's surface need not be reduced to the geoid (which differs from the 
surface of the ocean by the "sea surface topography, " on the order of 1 m), and 
secondly, satellite altimetry as additional observed data fits very neatly into the 
process. 

For reasons to be explained later, the analysis here is essentially devoted 
to investigating the effects on the undulation due to erroneous potential coefficients, 
as well as the lack of higher-degree coefficients. The discussion is further lim- 
ited to the case in which the earth is approximated by a sphere. Moritz (1974) 
treats the effect of the earth's ellipticity; Christodoulidis (1976) thoroughly 
studied the corrections for the atmosphere which are applicable also in this case; 
and finally, the terrain correction is negligibly small on the oceans. 

2. The Concept in Principle 

The known orbit of a satellite determines its altitude above the reference 
ellipsoid. The altimeter measurement a establishes the sea surface height 
h above the ellipsoid (see Figure 1). This height is the sum of the geoid undu- 
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lation N and the sea surface topography t , since it is generally assumed that 
the three surfaces depicted in the figure are parallel locally. 



Gravity g is measured on the sea surface, e.g. at point P, while the normal 
gravity y P at P is computed with the aid of the sea surface height: 

y, = y, + h (i) 

where the normal gravity y q on the ellipsoid and its radial derivative are known 
quantities. Thus, the sea surface height as obtained from the altimeter measure- 
ment serves the dual purpose of providing normal gravity at P , as well as the 
sea surface topography for a known undulation. 

We note that the gravity measurements on the surface could be reduced to 
the geoid if the heights of the sea surface above the geoid were known; but pre- 
cisely these heights are to be determined. Therefore, in practice it is simpler 
conceptually to obtain gravity disturbances (aided by altimetry data) on the ocean 
surface, than gravity anomalies on the geoid, the latter being required for any 
use of Stokes ' formula. 


From the definition of the gravity disturbance, 6g = g P - Tp , and the 
disturbing potential T, there is the basic relation (Heiskanen and Moritz 1967, 
P* 85): 

a g = - - < 2 > 

From this and with a spherical approximation, Moritz (1974) derives 


T = — 
4tt 


5g S ( 0) dcr 


where cr is the unit sphere, and 6g is given everywhere on the earth's surface 
(assumed to be a sphere of radius R). S ( 0) is an analogue to Stokes' function, 
defined by 
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(4) 


- , V 1 2n + 1 _ , 0 „ . 0 

S(<&) = L nTl Pn (cos ^ = csc ~2 ( 1+CSC Y ) 

n =0 

The argument lj) is the usual central angle between the point at which T is 
computed (on the sea surface) and the surface element dcr. The closed ex- 
pression for the infinite sum is derived by Hotine (1969, p. 311). Note that the 
sum in (4) includes zero- and first-degree harmonics; whereas, these are 
absent in Stokes' function. 

Let y be the normal gravity in the ellipsoid. The difference between the 
disturbing potential on the geoid, T = N 7 ( Bruns' formula), and the disturbing 

potential on the sea surface (equ. (3)) is approximately t 6 g and is neglected 

(for t=2m, Sg = 40mgal, N = 30 m, it is only .003% of T). Hence, the 
geoid undulation in terms of gravity disturbances is 

N = -JL- f f 6g S(0) da (5) 

4 rry J J 
a 

It is mentioned again that this, as Stokes' integral, is only a spherical approx- 
imation, so that y becomes an average value of normal gravity. 

The integral in (5) is unlike the usual Stokes' integral, however, in that it 
incorporates the zero-degree and first -degree potential harmonics. Whereas, 
the gravity anomalies in Stokes' integral must average to zero (so that N is the 
geoid undulation with respect to the given reference ellipsoid ) and have no first- 
degree term in their harmonic expansion, no such provisions are necessary for 
the gravity disturbances in (3) or ( 5) (see also Hotine 1969, p. 317). Never- 
theless, since the error analysis will require a decomposition of the undulation 
N into spherical harmonics, the zero-degree and first-degree contributions will 
be singled out. Certainly, if the global gravity disturbance data are known to 
possess neither zero- nor first-degree harmonics, then S($) could be modified 
by also excluding its first two harmonics ( see the end of section 3 and Appendix C). 


3. Truncation Coefficients 

To compute the undulation from an integral equation such as (5), gravity 
must be known, in theory, at every point on the earth. In practice, of course, 
such detailed gravity information is not attainable. Instead, only a finite set of 
measurements is available in the vicinity of the point of computation. Limiting 
the integration in (5) to this area is not sufficient, however, since geoid undula- 
tions are predominantly the consequence of the global features of the gravity 
field (as, for example, an inspection of the degree variances of the undulation 
covariance function shows). Therefore, the influence of the "remote zones" 
must be taken into account, and this is usually accomplished by separating the 
contributions of the inner and outer zones as follows: 
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■J j6gS(0) 


dcr + 


J j6gS(0)< 


cr-CTc 


where cr c is a spherical cap with radius 0 O , cr-cr c being the "cap'* with 
radius rr- 0 O . The evaluation of N 2 is predicated on a given set of potential 
harmonic coefficients, and therefore, should be rewritten in the form of an 
infinite series. A way to proceed follows the lines of development given by 
Heiskanen and Moritz (1967), sec. 7.4, starting with the definition 


S (llM = f 0 ’ i£ 0< * < ^o 

1<V) IS (ill), if 


S x ( 0) is piecewise continuous in the interval 0 < 0 < it , and hence, this 
function is expandable in a series of Legendre polynomials: 

CO 

s l( 0) = y k n Q n P n (cos 0) (9) 

nl= 0 

The coefficients k n Q n are determined by multiplying both sides of equation (9) 
by P n (cos 0) and integrating over 0 < 0 < tt (taking advantage of the orthog- 
onality of P n ) : 

k n Qn = J S l( 0) P n (cos 0) sin 0d 0 (10) 


With the choice 


k„ = 


2n + l 


equation ( 10 ) becomes, upon substituting ( 8 ), 

r n - 

Qt,(0o) = S(0) P n (cos 0) Sin0d0 (11) 

J 0o 

Let coordinates (0,a) be spherical coordinates on a sphere where the pole is 
the point ( 9 ,A) at which N is to be computed. Then, from ( 6) and (7), 

Na = 7“ rj J. 6g ( 0 , a ) §■( 0) Sin 0 d0 d a: (12) 

417 y J o 

which in view of ( 8 ) and ( 9 ) is 

U V 2n-i-l _ p 2 tt p rr 

N 3 - \ n ^g(0,«) P n (eos0)sin0d0da (13) 

4lT ' n =0 * J 0 ‘'O 

The gravity disturbance can also be expanded in a, series of Laplace harmonics 
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(let 0, X be the spherical coordinates of the point at which N a is computed): 

CO 

6g(0,X) = Y 6g n (0,X) (14) 

n — 0 

These Laplace harmonics, conversely, can be found by integrating the gravity 
distrubances (let (0,X) be the pole of the earth and use equ. ( 1-71) of Heiskanen 
and Moritz 1967 ) : 

2 n+ 1 p 2 tt pTT 

6gn ( M) = J 0 j Q 6g ((i>,a)P n (cos(i)sin0d0da (15) 


A comparison of equations (15) and (13) shows that 


CO 



(16) 


To achieve the objective of expanding N s in a series of potential coefficients, 
we write 


n = 0 
n 

with T n = [ C nm cos mX + S nB sin mA ] P nn (cos0) 

P nn (cos 0) is a fully normalized Legendre function; C na , S nn are fully 
normalized potential coefficients, with C nn being the coefficients after the 
normal gravity field has been removed; and kM is the product of the gravita- 
tional constant and the mass of the earth. Finally, with (Moritz 1974, p. 48) 


• (17) 


6 Sn = T n and y = ^ (18) 

K R 

equation (16) becomes 

co n 

N 2 (0,X) = ~y (n+l)Q n y [C na cosmX + S na sinmX] P nB (cos 0) (19) 

n = 0 

Remembering that Q n is a function of the cap radius $ 0 , it is noted that for 
$ o = 0, N a = N . For 0 O > 0 , the coefficients Q n , which are analogous to 
Molodenskii's truncation coefficients ( Molodenskii, et al. 1962, p. 147), 
actually effect an increase in the rate of convergence of this series, as the 
contribution of the finer structure of the gravity field (the high-degree harmonics) 
is essentially accounted for by N x in equation (7). Of course, this increase in 
the convergence rate depends on the size of $ 0 . 

The difference between Stokes' integral and its analogue in terms of 
gravity disturbances ( equation (5) ) has already been elucidated in the previous 
section. This difference has a direct bearing on the truncation coefficients Q n . 


It is again noted that as opposed to Stokes’ function the series expansion of 
S (0) ( equation (4) ) includes zero- and first-degree terms. If 6g 0 and 5g x in 
the series expansion of <5g (equation (14)) were zero, then the integral formula 
(5) in concept would not change; bat in order to be consistent with Stokes' in- 
tegral, the first two harmonics of S(0) could be excluded. This modification 
of the kernel function evidently induces a change in the truncation coefficients 
Q n ( see also section 10 ). Since the formulation of the truncation coefficients 
as presented above yields the higher error estimates in the analyses of section 
7 ( see also DeWitle 1966), the alternative approach using the modified kernel 
function S* (0) = S { 0) - 1 - -|-cos 0 is pursued only in Appendix C. 

4. The Evaluation of Q n 

There exists a handful of methods to evaluate the coefficients Q n ; not all 
are equally accurate and efficient- Direct integration by analytic means seems 
plausible for small n (as Molodenskii has done for his coefficients). A recourse 
to numerical integration for larger n is unsatisfactory in view of the limited 
attainable accuracy and the singularity of S(0) at 0 =0 (see equation (4)). It 
is possible also to determine a finite sum for each Q n since P n is expressible 
in this way. But the volume of computations (coupled with a loss of accuracy) 
for large n is prohibitive. A desirable formula for Q n would be in the form 
of a recurrence relationship, such as Paul (1973) has derived for Molodenskii's 
coefficients by utilizing the series expansion of Stokes’ function (Hagiwara (1976) 
also formulated a very elegant recursive solution). 

In the case of the coefficients Q n , the integral (11) actually admits to a 
simplifying integraJon by parts, and by borrowing several results from Hagi- 
wara (1976), a relatively simple recurrence formula can be derived. The appen- 
dix gives a full account of this derivation for ip 0 ^ 0. This solution of Q n seems 
to yield quite accurate and very efficient evaluations even for small values of 
0 O . Although S(ipo) has a singularity at ip o =0, Q n is nevertheless well- 
defined there; for indeed, a comparison of (9) and (4) shows (noting that when 
0 o = O, S(>p) = S x (0) for O<0< tt ) 

Qn (0) = • , n a 0 (20) 


5. The Global Error Estimates 


The total geoid undulation at a point can be expressed as an infinite series 
of spherical harmonics by combining Bruns' formula (N = T/y), and equations 
(17) and (18): 


N 



6 gu 
n + 1 


( 21 ) 


Considering equations (7) and (16), (21) may be divided as follows; 
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N 

Ni 

Na 


N x + Na 

ft, (nTl - HO 6 s» 


_R 

2 y 


[Qo 5g 0 + Qi 6gJ 


IV. 


27 L 

n =2 


Qn Sgn 




> (22) 


/ 


It is usual to assume in gravimetric determinations of the geoid that the centers 
of mass of the reference ellipsoid and the earth coincide. Then 6g a = 0. How- 
ever, their total masses may not be equal; in this case, the zero-degree term 
does not vanish and from (18) and Heiskanen and Moritz (1967, p. 98) 


R - . _ _R_ — To _ Qo k 6 m 

2 y Qo °go - 27 Qo R “ 2 7 R 


(23) 


where 6 m is the difference in the total masses. N x also contains a zero-degree 
term, but it is not necessary to treat it separately if Nx is evaluated according 
to the integral expression of equation (6). To be observed is that the difference, 
6W, between the ellipsoidal and geoidal potentials has so far been ignored. If 
6w ^ 0, then by the generalized Bruns’ formula ( N = ( T - 6W )/y) the total con- 
stant part of N is k6M/‘yR - 6w/y. 

The first infinite sum in (22) is the contribution to N from the gravity 
field within the cap cr c (see equ. (6) ). To obtain a high degree of accuracy in 
N i , one usually visualizes a sufficiently dense set of measured gravity and 
altimetry values in cr c which are then processed in some systematic fashion 
resulting in an approximation to the first integral of (6). The errors in this 
approximation arise from errors in the gravity and altimetry data itself (pro- 
pagated error), as well as from an unavailability of continuous data over the 
cap (discretion error). Christodoulidis (1976) investigated two methods which 
could provide estimates for the errors in N x when it is formulated (analogous 
to equ. (6) ) as an integral of gravity anomalies . He concluded that these meth- 
ods are computationally unfeasible and proceeded to estimate the corresponding 
errors in the so-called ’’frequency domain’’, that is, when N x is expressed in 
a spectral representation as in equation (22). The structural differences between 
the formulations of the geoid undulation in terms of gravity anomalies and in 
terms of gravity disturbances are practically nonexistent; and hence, it will 
be assumed that in the present context, the same type of numerical difficulties 
would occur in the error estimation of N x . Moreover, and for the same rea- 
son, it is obvious that a similar procedure as developed by Christodoulidis can 
lead to error estimates for N x in the frequency domain. Relying on the argu- 
ment that the local behaviors of gravity anomalies and gravity disturbances are 
not inherently dissimilar, the derivations and evaluations of these error esti- 
mates are omitted with the assumption that no significant deviations from Chris- 
todoulidis' numerical results would be expected. 



To put some rigor into the underlying argument, we digress briefly and 
derive the covariance function for gravity disturbances. For the disturbing po- 
tential, the covariance function is given by (Moritz 1972) 


K(P,Q) 



P n (COS 0 ) 


(24) 


where 0 is the spherical distance between the points P and Q , r P , r Q are 
the radius vectors to P and Q , and the c n are the degree variances for gra- 
vity anomalies. By the law of propagation of covariances (Moritz 1972) as applied 
to the relation (2), the covariance function for gravity disturbances is 


Br3 “1 v/n+1 

D < p - Q > = -57 P t^ K(P ’ Q) j = l \—- 1 


R 


2 . n+ 2 


n=2 


r P r Q 


■) Bi(cos0) (25) 


By inspection, the degree variances for gravity disturbances are 

,3 

2 (26) 


fn + 1 \° 

dn ~ Vn - 1/ ° n 9 n ~ 


For large n , that is, for high frequencies, the difference between the degree 
variances d n and c„ is negligible, which supports the contention that the local 
variations of gravity anomalies and gravity disturbances are quite similar. 

Having made these comments, the remaining part of the report will be 
concerned only with errors in the contribution to N from the gravity field out- 
side the spherical cap c r c . Furthermore, only differences in undulations are 
treated, these being used in the definition of a relative geoid. This also effec- 
tively eliminates the need to consider errors in the zero-degree (constant) term 
of the undulation. Since N a (see equations (6) and (22) ) is to be computed from 
a given set of potential harmonic coefficients, it is natural to estimate the errors 
in the frequency domain; that is, we start with N a as defined by equation (22). 

Let P and Q be two points on the earth (sphere) at which Na is deter- 
mined. The error incurred in the difference AN a = N 2 P - N 2 q is 

Ac = ( Na P - Naq ) - ( Nap ~ Naq ) true (27) 

or Ac = e p - e Q (28) 

where e P and €q are the errors in N 2 P and Naq . The error e P is composed 
of two (presumable uncorrelated) errors; namely, the error e lP due to erroneous 
potential coefficients (commission error) and the error e ap resulting from the 
lack of higher-degree coefficients (omission or truncation error): 

Ae = Ae x + Ae a = ( e 1P - e lQ ) + ( e aP - e 2Q ) (29) 
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The error estimates that are studied in this report are global averages of 
the errors, the average being defined by the operator M: 


M(* ) 


— r r 

4tt J v 
a 


(•) da 


(30) 


This so-called ''space average" (Moritz 1978) serves also to define the covariance 
functions, such as equations (24) and (25). 


Consider first the error e 2P induced by the unavailability of harmonic co- 
efficients of degree greater than n ; from (22) it follows that 

OO 

e 2P = -Y Q n <5g n (31) 

Z / o _ 
n = n +1 

This expression is quite analogous to the expansion of the gravity disturbance into 
spherical harmonics as in equation (14). Recalling the covariance function for 
6g, equations (25) and (26), it is immediately evident that the covariance func- 
tion for C s (that is, the global average of the products e 2P e 2Q for points P, Q 
which are separated by a fixed spherical distance i/j Pq ) on the earth is given by 


M(e 2 p£aq) “ 


£ Qn 2 d n P n (cos ip PQ ) 


(32) 


n = n + 1 


When P and Q coincide (^> Pq = 0), then M(e 2 2 ) is the variance of the error 
e 2 (independent of the location of point P on the earth) 


M(€ 2 ) 



(33) 


On the other hand, the variance (global mean square) of the difference error 
Ae 2 is 

M((A<f 2 ) 2 ) = M( e| P ) + M( c 2 q) - 2M( e 2P e 3Q ) (34) 

which when equations (32), (33), and (26) are substituted becomes 

CO 

rn /tt+lv 

M((A€ 2 ) 2 ) = C n (! “ Pn(COS ^pq) ) (35) 

' n = n + l 

The degree variances c n in this equation refer to the mean earth sphere of 
radius R . If, on the other hand, they are computed according to the model 
of Tscherning and Rapp (1974), for example, then the sphere of reference is 
the Bjerhammar sphere (the sphere entirely embedded within the earth). To 
rectify equation (35), the factors a n + 2 = (Re 3 /R 3 ) n+2 are introduced 
(Jekeli 1978). Henceforth, unless otherwise stated c n is understood to refer 
to the sphere of radius Re ; then 


■ 9 " 



( 36 ) 


2 00 2 

mi = M((Ae a ) a )=^^ ^" 3 £rf) e » ^” + 2 (l-P»(° os <M> 

ma is the global estimate of the error in AN 2 due to a lack of higher-degree 
potential coefficients. This equation is obviously valid for any value of R> R B , 
but we restrict ourselves to the earth's surface, where R = 6371 km. 


For comparison, the corresponding error estimate derived by Christodou- 
lidis (1976) is ^ 

(m ' s f = Q n 2 c n CT n+2 (l-P R (cos </) Pq )) (37) 

where the Q n ( , n = 2, 3, ... . are Molodenskii's truncation coefficients. 


The remaining error e i , caused by noise in the available potential coeffi- 
cients can be estimated similarly. Let C nm , S nm be the true, fully normalized 
harmonic coefficients (see also equ. (17) )and let C„' B , be the corresponding 
coefficients that are available from combined, global gravity and satellite data 
(see also section 6). Then explicitly, using equation (19), the error in N 2 can 
be expanded in the series (note that the zero-degree term has been omitted) 


e-i = 


OO n 

“ |^(n+l)Q n ^[(Ci B -C nm )cos rnX + (S; m -S nra )sinmX]P nn (cose) (38) 


n=2 


m ~ 0 


If the covariance function on the earth of this error component is written in the 
form 


M ( € 1P c lQ ) 


CO 

Y Pn (COS 0 PQ ) 
n=2 


(39) 


where tf) PQ is the spherical distance between points P and Q , then the degree 
variances £ n are given by (seee.g. Moritz 1972). 


2 n 

= ^_(n + l) 2 Q n 2 V [(C' B -C nffl ) 2 +(S; ffl -S nm ) 3 ] , n^> 2 (40) 

m = 0 


Suppose that n is the maximum degree to which these coefficients have been 
determined, then obviously we can take 


c = 0 = C nn , S n T m = 0 = S nm , for n > n ( 41 ) 

since the error due to neglected coefficients past degree "n has already been 
accounted for. 

The true values C nm and are unknown, therefore, the differences 
appearing in equation (40) must be replaced by the corresponding estimated 
standard errors of the coefficients, it is shown later that under certain conditions, 
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we can assume equal standard errors for C n T ffl and S n ! j, for each degree n : 

mc na = m B = ms na , 2 < n ^ n (42) 

A 

Then if £ n denotes the estimate of £ n based on (42), 

£» = (n + l) 2 (2n + l)m n 2 Q n 2 , 2 n ^ n (43) 

A 

and in view of (41), £ n = 0 for n > n ; which yields finally the following 
estimate of the covariance function (£9) on the earth: 

2 \ 

M(€i p ^ ) = “.^(n + l) 2 (2n+l)m n 2 Q n 2 P n (cos(iip Q ) (44) 

n=2 


As in the case of Ae 2 (see equation (34) ), the estimated global mean square 
of the difference error Ae x is 

MUA^) 3 ] = 2[M(e 2 3 ) - M(e 2p e 2Q )] 

n + l) 2 (2n + l)m n 2 Q 8 (l-P n (cos(()p q )) (45) 



Again, if the potential coefficients refer to a sphere of radius Ra ( in the sense 
that R in equation (17) is replaced by Re ), while the undulation differences 
are determined on a sphere of radius R , consistency is restored by multiply- 
ing 4 n by cr" ( see Appendix B ) where cr = R b 3 /R 3 : 


m. 


2 . 

- M^Ae,) 2 ] = ^-^Q n 3 (n+l) 2 (2n + l)m n 2 cr D (l-P n (cos^p q )) (46) 


mi is then a global estimate of the error in A Ns that results from the random 
errors in the given potential harmonic coefficients. Christodou'lidis (1976) 
derives the following estimate of the corresponding error in terms of the 
Molodenskii truncation coefficients Q n , written here for comparison and later 
reference: 

2 n 

(mi) 3 = ~^Qn 3 (n-l) 3 (2n+l)m 3 a n (l-P n (cos^ Q )) (47) 

n =2 


As briefly remarked earlier, the estimates m x (equation (46)) and m 2 
(equation (36)) may be regarded as uncorrelated errors; and therefore, the 
total estimated error in AN 2 is 

m = S mj 3 + m a 2 (48) 

Some final notes with respect to the computation of these errors are 
appropriate. Both the commission error m-L and the truncation error m 2 
are subject to further approximation in the type of analysis (simulation study) 
to be conducted here. The degree variances c n for gravity anomalies and 
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m n 2 for the errors in potential coefficients are obtained through simple models. 
The ( global) model for c n is usually determined on the basis of "observed” 
degree variances of low degree (e.g. Tscheming and Rapp 1974). However, 
the sum in (36) requires degree variances for the high frequencies of the gravity 
field; this makes the choice of the model somewhat problematic in that it should 
actually reflect primarily these high-degree variations of gravity anomalies. 

The model for m n 2 , on the other hand, depends to some extent on the method 
by which the potential coefficients are determined. Christodoulidis (1976) 
constructed a model first tailored specifically to agree with the accuracies of 
the GEM 6 solution ( (16, 16) combined solution) and then generalized in a 
natural way for similar higher-degree solutions. We shall see later, in the 
light of newer and improved data, that this model may no longer be entirely 
appropriate. One further approximation in the evaluation of m^ 2 replaces the 
upper limit 05 of the sum by a number, but with a fast computer there is no 
difficulty in making this number sufficiently large. 


6. The Models for c n and m n 3 


Little need be said in regard to the c n -model - only that the model of 
Tscheming and Rapp (1974) was used here exclusively. The effects of a different 
model (the one proposed by Moritz (1976)) on the truncation error will be inves- 
tigated briefly in section 9. 

The model for m n 2 that was selected by Christodoulidis (1976) has the 
form ( F is a constant ) 





y 3 ( n - 1 ) 2 ( 2n-i-l ) 


2 s n ^ n , n > 16 (49) 


Noting that the degree variances for gravity anomalies are given by 

n 

c n - y 2 (n-l) 2 £ (C n 2 + S 2 ), n s 2 (50) 

in =0 

it is readily seen that when n = n in (49), F 2 is the relative error of the 
potential coefficients of degree n . Therefore, this model was designed on the 
basis that every solution for the potential coefficients yields approximately the 
same relative error at the maximum degree (Christodoulidis chose F = .8 , 
implying a 64% relative error at If). 

Potential coefficients can naturally be determined from purely terrestrial 
gravity data, as well ( Rapp 1969 ) ; 



1 

4n y (n- 1) 



(51) 


Given a global distribution of mean gravity anomalies, the integral must be 
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approximated by some type of sum, e. g. 



4tt y ( n - 1 ) 


y 


_ f cos mX 

4gP„.(cos6){ sinml 


Aa 


(52) 


where Ag is the_ observed mean gravity anomaly for the block Acr that is 
centered at ( 0 ,X ) . The potential coefficients are computed through a least 
squares adjustment. It is not difficult to show (Rapp 1969) that the standard 
error ( i. e. the square root of the diagonal element of the covariance matrix) 
of the coefficients of degree n is given by 


m 


n 



A0 AX 

V-i) 


7T 


n :> 2, A 0, AX in radians 


(53) 


This formula is derived with the assumption that the mean gravity anomalies 
are observed over the whole earth with the same accuracy . Also, the 
error due to the fact that these observations constitute a noncontinous data 
set is ignored. It will be assumed that any such fault in the standard error 
above is incorporated into the overall uncertainty of the modeled error variance. 


■ For a given solution of the coefficients C„ B , S nm all quantities (except n) 
on the right side of (53) are constant; thus, we may say that the error in the 
potential coefficients of degree n > 2, as determined from a global set of mean 
anomalies, is inversely proportional to (n-1). Hence, the following simple 
model suggests itself for this case: 


m n 


K 


n - 1 


K = constant , 2 £ n £ n 


(54) 


Table 1 below compares die models given by equation (49) and (54) against 
the errors m n of the solution obtained from GEM 9 plus 1° x 1° mean anomalies 
(Rapp 197S, private communication) up to n = 12 . In equation (49), we take 
F = . 8 , n = 20 ( GEM 9 is a ( 20, 20 ) satellite solution), 


k 

Table 1. Comparison of Models (49) and (54) with m n 
(obtained for the solution using GEM 9 plus l°x 1° mean anomalies) 


n 

m n x 10 6 

m n x 10 fa , model (49) 

m n x 10 s , model (54) 

2 

.0027 

.0080 

.0581 

3 

.0066 

.0155 

.0291 

4 

.0046 

.0138 

.0194 

5 

.0072 

.0136 

.0145 

6 

.0057 

.0137 

.0116 

7 

.0067 

.0139 

.0097 

8 

.0057 

.0142 

.0083 

9 

.0057 

.0145 

.0073 

10 

.0051 

.0148 

.0065 

11 

.0050 

.0151 

.0058 

12 

.0044 

.0154 

.0053 
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and the c n -model of Tscherning and Rapp (1974) (see equation (59)). In 
equation (54), K = .0581 x 10~°. This value of K is the average of the values 
of IC implied by the standard errors from n = 13 to n = 60 of potential 
coefficients computed from 1° x 1° gravity data alone (Rapp 1978, private 
communication). Table 2 shows the relevant information by comparing the 
resulting model (54) with the standard errors of the potential coefficients. 


Table 2. Comparison of model (54) ( K = .0581 x 10 -s ) with m n obtained for 
the solution using 1° x 1° mean anomalies ( Rapp 1978). 


n 

m n x 10 s 

m n x 10 s , model (54) 

n 

A c 

m n x 10 6 

m n x 10 s , model (54) 

13 

.00483 

.00484 

37 

.00162 

.00161 

14 

. 00446 

.00447 

38 

.00157 

.00157 

15 

.00414 

.00415 

39 

.00153 

.00153 

16 

.00386 

.00387 

40 

.00149 

.00149 

17 

.00362 

.00363 

41 

.00145 

.00145 

IS 

. 00341 

.00342 

42 

.00142 

.00142 

19 

.00322 

.00323 

43 

.00138 

.00138 

20 

.00305 

.00306 

44 

.00135 

.00135 

21 

.00290 

.00290 

45 

.00132 

.00132 

22 

.00276 

.00277 

46 

.00129 

.00129 

23 

.00264 

.00264 

47 

.00126 

.00126 

24 

.00252 

.00253 

48 

.00124 

.00124 

25 

.00242 

.00242 

49 

.00121 

.00121 

26 

.00232 

.00232 

50 

.00119 

.00119 

27 

.00223 

.00223 

51 

.00116 

.00116 

28 

.00215 

.00215 

52 

.00114 

.00114 

29 

.00207 

.00207 

53 

.00112 

.00112 

30 

.00200 

.00200 

54 

.00110 

.00110 

31 

.00193 

.00194 

55 

.00108 

.00108 

32 

.00187 

.00187 

56 

.00106 

.00106 

33 

.00181 

.00182 

57 

.00104 

.00104 

34 

.00176 

.00176 

58 

.00102 

.00102 

35 

.00171 

.00171 

59 

.00101 

.00100 

36 

.00166 

.00166 

60 

.00099 

.00098 


The GEM 9 solution combined with the 1° x 1° terrestrial data yields one 
of the most recent and realistic accuracy estimates for the potential coefficients. 
This solution supplies coefficients up to a maximum degree n = 60 ; but 
solutions up to n = 180 have been done (Rapp 1979, private communication). An 
inspection of Table 1 reveals the inadequacy of model (49). Different values for 
the model parameters F and n could conceivably diminish the discrepancies, 
but such changes in the parameters are not easily justified from a theoretical 
standpoint since the model is purely an empirical construction. Consequently, 
the model (49) is not studied further in this report. Table 2 exhibits the ex- 
cellent agreement between the model (54) and the errors of potential coefficients 
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( for 13 £ n ^ 60 ) that are derived from terrestrial data. However, model 
(54) compares not as well with the low-degree ( n ^ 6 ) errors obtained through 
the combined solution ( see Table 1) which are favorably influenced by the "a 
priori” coefficient errors from the satellite solution. In specialized investiga- 
tions, it is not unfeasible to use a model such as (54) for the high-degree errors 
and the actual standard deviation, m n , as given in Table 1 (first column) for 
the errors at the low-degree end of the spectrum: 


m 


n 


r m n 

{-Sr 

n- 1 


2 ^ n ^ n 2 
n x < n £ h 


(55) 


7. Computational Analyses of the Errors 

The radical dissimilarity between model (54) and model (49) which was 
used by Christodoulidis (1976) warrants some new tests here with his error 
equations (37) and (47). (It should be remarked that these equations are com- 
pletely valid provided that gravity anomaly data referred to the geoid are avail- 
able within the cap a c .) Since m x and m 2 are the square roots of squared 
errors averaged over the earth, they are designated root mean square ( RMS ) 
errors in the figures. They are also called simply error estimates in the 
following discussion. To avoid unnecessarily cumbersome language, the 
following designations are introduced: 

Method A = the method of computing error estimates of 
AN based on gravity anomaly data within the 
cap a c (equations (37) and (47) ) 

(56) 

Method B = the method of computing error estimates of 
AN based on gravity disturbance data within 
the cap a c (equations (36) and (46) ) 

The error estimates of AN a may be investigated with respect to variations in 
the separation (f) PQ , in the cap radius ij) 0 , in the maximum degree n of the 
potential coefficients, and in the model (54) itself. 

Method A was used to determine the general characteristics of the AN a - 
errors as they are implied by the model (55) with K = .0581 x 10 ~ e , n 3 - 12, 
and m n given in Table 1. Figure 2 depicts, for = 30°, the commission and 
truncation errors in AN a as functions of the separation i/fpq and the maximum 
degree n. For these and all subsequent evaluations of the truncation error, 
the upper limit >» in (36) and (37) is "approximated” by 3000 at which point a 
precision of at least .1mm seems to be attained. Evidently, since the trunca- 
tion error is independent of the standard errors of the available potential coeffi- 
cients, the corresponding curves trace the same values given by Christodoulidis 
(1976, Table 7, n = 30). The commission errors, however, increase with 
increasing n, as opposed to the corresponding decrease in the commission 
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Figure 2. Commission (Solid line) and Truncation (Dotted line) 
Undulation Difference Errors Using Current Potential 
Coefficient Error Estimates and Assuming Perfect 
Gravity Anomaly Data Within a 30° Cap. 
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errors of Christodoulidis ( see his Table 6). These contrasting characteristics 
are obviously the direct consequence of the difference in the error models, 
specifically, because the model (54) or (55) does not stipulate the same relative 
accuracy at all values of n. As seen in Figure 2, the commission errors do 
appear to be bounded ( as n varies ), at least for larger values of $ PQ . It is 
observed that the truncation error loses its dependence on the separation 0 Pq 
as the number of available potential coefficients increases, Contrarily, the 
commission error varies significantly (but not without bound ) , regardless of 
the value of n , as the spherical distance between points P and Q increases. 

The root sum of squared truncation and commission errors of Figure 2 is 
presented in Figure 3 for both Methods A and B. Figures 4 and 5 are analogous 
to Figures 2 and 3 in that the same error model (55) is used in the computation 
of the errors. But here, the errors are shown as functions of the separation 
0 p q and of the cap radius ipo, n being fixed at 180. A ’’cap radius" of (£> o = 0 
signifies that there is no cap and implies that the undulations are determined 
from potential coefficients only. Again, the truncation errors are practically 
independent of the separation 0 PQ (since n = 180, cf. Fig. 2). The commission 
errors exhibit a similar type of dependence on 0p Q as in Figure 2; and they 
must plainly increase strongly as the size of the cap a c decreases implying 
less and less detailed gravity data near the computation points. The total ( RMS ) 
error as a function of 0 PQ and tj) 0 , shown in Figure 5, in this case resembles 
the dominating commission error. 

It is noted that for ifi 0 = 0 (no cap), the errors according to Method B are 
equal to those of Method A since no terrestrial gravity measurements "within 
the cap" enter into the computations of the undulations. In fact, the error equa- 
tions (36), (46) are identical to (37), (47) when Q n (0) = 2/(n+l) and Q n (0) = 
2/(n-l) are substituted (see also Figure 7). 

Having established the general characteristics of the commission and 
truncation errors in Figures 2 and 4, the remaining investigations, concentrate 
on the total (RMS) errors in AN 2 with n fixed at 180. In addition, the potential 
coefficient error model (54) is now used for the entire spectrum, n = 2, . . . , 180. 
More explicitly, this model can be written (with y - 979800 mgal, A0 = AX = 

1° = .01745 rad. in equation (53) ) as 

0.00502 x 10~ 6 

m n = e. (57) 

n - 1 As * v 1 

where e^ s as a new variable in the analysis is the average error in the mean 
1° x 1° gravity anomalies. Figure 6 displays two sets of errors computed by 
Method A, for 0 o = O° and 0 o =3O°. The curve for $ o = 0°, = 1 mgal 

essentially traces the truncation error which in this case almost completely 
overshadows the commission error; and therefore, it represents the maximum 
attainable accuracy in undulation differences which are derived solely from 
potential coefficients up to degree 180. Recall that the standard errors as 
modeled by (54) for small degree are much larger than what might be expected 
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Figure 3. Total RMS Errors for Geoid Undulation Differences 

Assuming Perfect Gravity Data Within a 30° Cap Using 
Gravity Anomalies (Method A) and Gravity Disturbances, 
(Method B), and Current Potential Coefficient Error 
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Figure 6. 
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from a combined solution ( see Table 1, columns 1 and 3 ). We find that this is 
the reason for the steeper increase of the ( = 0, = 10 mgal ) - curve at 

the larger values of t/) F , Q when compared to the corresponding curve of Figure 5 
(where K = .0581 x 10 -8 implies = 11.6 mgal). Finally, Figure 6 includes 
curves for errors in AN 2 when the undulations are referred to a perfect 30- 
degree potential reference field. These errors are computed by replacing the 
lower limit of the sum in equation (37) by n = 31. Results ( not shown) similar 
to those of Figure 6 can easily be produced also for Method B ; the same general 
characteristics should be expected. 

Figures 3 and 5 have shown that the errors, whether through Method A 
or through Method B, possess a parallel dependence on both the maximum degree 
n and the separation tf)pq ; Figure 7 finally compares Methods A and B in regard 
to the dependence on the cap radius • These curves represent the errors only 
for a particular separation, ^pq = 1500 km. Other values of 0 Pq would not alter 
the basic interrelationships shown here. 

As observed in Figure 3, the total errors of Method B are greater than the 
corresponding total errors of Method A when the cap radius is 30° and the model 
(55) is chosen. However, Figure 7 (where the model (54) is used) and Figure 5 
both show that this is not normally the case. 

It should be strongly emphasized that, except for the case ip 0 ~ 0 , all 
figures depict estimated, average errors associated with only the outer zone 
contribution, AN a , to the total undulation difference. Recall that iJj 0 equal to 
zero implies that AN 2 = AN is the total undulation difference. Up to this point, 
no numerical description has been given to the errors that are encountered within 
the spherical cap. 

If the data of this inner zone were continuous and errorless, thus yielding 
exact values of AN 1S the error in AN 2 , of course, would also be the error in 
the total undulation difference. In the absence of this ideal situation, two types 
of errors in the gravity data enter into the computation of N x ; namely, the 
discretion error resulting from the lack of continuous data, and a propagated 
error due to imperfect measurements. The latter includes the influence of 
errors in altimeter measurements ( see section 2). 

It is important to realize that the principle of truncation functions as 
utilized here, and in general, is based on the assumption that the gravity ma- 
terial in the vicinity of the computation point ( i.e. within the spherical cap) is 
more dense and accurate than in the remote zones. Indeed, for a global, uni- 
formly accurate field of mean anomalies it is pointless, if not invalid, to com- 
pute first potential coefficients ( and subsequently N 2 ) from this global set and 
then N-l from the same mean gravity anomalies in a spherical cap. In this 
special case of uniform coverage with no possibility of local densification, we 
must choose 0 O = 0 „ 



Figure 7. Total RMS Errors for Geoicl Undulation Differences at 1500 
km Assuming Perfect Gravity Data Within a Cap of Radius 
fjio with Potential Coefficients Estimated to Degree 180 Based 
on Anomalies Given to an Accuracy of €^ g . 
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Rapp (1973) and Christodoulidis (1976) have examined these errors with 
respect to the inner zone. Christodoulidis formulates the errors in the differ- 
ences AN X based on Method A. The argument of section 5, that the errors of 
AN X do not differ substantially in going from Method A to Method B, can be used 
to obtain average estimates of the discretion and propagated errors of An x by 
quoting results from Christodoulidis (1976). Only an example is given here to 
illustrate the orders of magnitude that may be expected. 

Assume that potential coefficients are available to degree 180 as obtained 
from a global field of 1° x 1° mean anomalies that have a uniform accuracy of 
±5 mgal. Further, suppose that two spherical caps, each having a radius of 
j/) 0 = 25° are covered by 30' x 30 T mean gravity disturbances with an average 
accuracy of ±1 mgal. Let's say that the cap centers are separated by 1500 km. 
The discretion and propagated errors associated with AN X ( see Tables 3 and 4 
of Christodoulidis, 1976) are 32.2 cm and 27.4 cm, respectively. From Figure 
7, the total error in AN 2 is 9.2 cm. Provided that AN X and AN 2 are uncor- 
related, the average, estimated error in AN, the total undulation difference, is 
43. 3 cm. 


8. Mean Undulation Difference Errors 


Thus far, the error analysis has been confined to differences between 
point undulations. The mean undulation function which represents an average of 
undulations over a spherical cap is characteristically derived from the point 
undulation function by introducing a smoothing factor /S n into its spherical har- 
monic expansion. Details of this procedure maybe found in Meissl (1971). Very 
briefly, the average undulation over a cap cr x of radius i}) 1 on the sphere can 
be expressed as 


N (&,\) 


1 

2tt ( .1 - cos ip x ) 



N( 9 , A) dcr 


(58) 


where 2 rr ( 1 - cos $ x ) is the area of the cap a x (on the unit sphere ) , and 
(’6, A) are spherical coordinates associated with this cap, for example, its cen- 
ter coordinates. The eigenfunctions of the integral operator above ( operating 
on N) are the spherical harmonics sin mX P nil ( cos 0) and cos mA P na ( cos0); 
while the corresponding eigenvalues are (for all orders 0 ^ m £ n, and n > 0) 

& “ sh < 59 > 


Therefore, when the spherical harmonic expansion of N ( equation (17) with 
N - T/y ) is substituted into (58), the mean undulation becomes (zero- and 
first-degree terms are omitted, see section 5) 


N(0,A) 


co n 

= R^0 n ^) (C,i„cos mA + S na sin mA) P nB ( cos 0) 

n= H a = 0 


(60) 
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The only difference between N and N is in the harmonic coefficients. This 
fact implies that we may write ( cf. equation (22) ): 


where 


N 


Ni + N a 

- fS(i£i - !«*)«. ««. 

n =3 

co 

N s = f-Y Q„ & 5g n 

^ ' n==3 



This aualogy can be extended to the point that the average, estimated 
errors in A N 3 ( = N aP - N aQ ) due to erroneous potential coefficients up to 
degree h and due to a lack of higher-degree coefficients are given by ( cf. 
equations (46) and (36), respectively ) 

a ” 

m 3 = + + (l-P*(cos tfV)) 

n — 2 

3 “ 3 

= 2^ l®" (^1 ) Sn 3 0„<J ” +2 ( 1 - P„ ( COS * P , ) ) 

_ kM 

y ~ IF 

It is assumed that if > 0, then ip } < ; that is, the mean undulation cap 

is smaller than the cap representing the inner zone. If ifi 0 ~ 0 r then the mean 
undulation is computed solely from a finite set of potential coefficients. Similar 
error equations for AN a can be derived on the supposition that the inner zone 
is covered with gravity anomaly data ( cf . equations (47) and (37), respectively): 


\ 


(62) 


J 


3 u 

( ml) 3 = -^-^Qn 3 (u-l) 3 (2n + l)^fm n 3 a n (1 -P n (cos </)*,)) (63) 


n = 2 


<m i) a « Q,f Af-c n a n + s (l-P n (cos 0 PQ )), y 

n='rT+l 


(04) 


Normally, the cap associated with the mean undulation is approximated by 
a block which has the same area. Therefore, fora 0 ° x Q x ° mean undulation, 
the corresponding cap radius \ii x is roughly (since 0, a v n $ x ) 

ip x = e x /s^r ( 65 ) 


where B x and ip x have tlie same units. For 0^ = 1°, tp x " 0?564; and for 
G x - 5°, ~ 2? 321. 

Mean undulations are presently envisaged to be derived from gravity anomaly 
data; and therefore, the estimated errors are displayed for Method A only. More- 
over, of the many combinations of the variables 0 O , n, , 0 X that could be in- 
vestigated, only a few of particular interest are shown in the figures. 
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For small values of 0 X , such as 1° or 2°, the attenuation of the smoothing 
factor jS n (equation (59) ) is rather slow for the low frequencies ( e.g. = 1°; 

,8 10 =.999, = . 8S3, /J 300 = . 240). Therefore, the effect of |8 n in the trun- 

cation error ( equation (64)) is much more pronounced than in the commission 
error (equation (63)). This is verified in Figure 8, where two pairs of curves 
(for = 1 mgal and e^s = 10 mgal) trace the estimated errors in the AN 2 
contribution to 1° x 1° mean undulation differences. Figure 8 should be com- 
pared to Figure 6 which depicts error curves for point undulation differences. 

If = 10 mgal ( and n = 180 ), the commission error dominates the truncation 
error and the total error generally falls by only 10% to 15% in going from point 
to 1° x 1° mean undulations. Note however, the more substantial decrease when 
the centers of the 1° x 1° blocks are close together ( 30% for 0pg = 100 km). 
Similar results for 9 X ~ 5° , e^ 6 = 10 mgal ( not shown) mostly show an additional 
decline in the errors of about 20% to 30% ( somewhat more for small 0 PQ ). 

The more dramatic effect is in mean undulation errors which reflect pre- 
dominantly the neglect of high-frequency gravity information. Evidence to this 
is given by the two curves of Figures 6 and 8 which correspond to the case 
0 o = 0, e = 1 mgal and show a difference of 50% between them. A further re- 
duction by 50% to estimated errors of 10 cm to 15 cm is obtained for 2° x 2° 
mean undulation differences. The curves in Figure 8 represent estimated errors 
in the total mean undulation difference AN for the case 0 O = 0 , but also when 
0 O = 20° if errorless differences AN X can be derived from the data of the inner 
zone. 


Figure 9 exhibits a comparison of mean undulation difference errors 
(0o = O°) with respect to various potential fields of reference. Little seems 
to be gained in accuracy when changing from a 20-degree to a 30-degree ref- 
erence field. Again, because the truncation error overshadows the commission 
error when = 1 mgal, the choice of the inference field, practically speaking, 
is immaterial in this case. 


9. The Effect of a Different c„- Model 


As mentioned previously, the estimated truncation error relies on the choice 
of the model for the anomaly degree variances. To examine this dependence, the 
errors with the model of Moritz (1976), 


« Oh 


n - 1 
n-s-A 


0 '- 


+ 3 


+ a 2 


n - 1 


(n-2)(n+B) 


cr 2 a+3 , 


n> 3 


( 66 ) 


arc compared to the errors with the model of Tscherning and Rapp (1974) which, 
so far, has been used without exception: 




n- 1 


( n-2)( n-HB) 


n > 3 


where ot a - 425.28 mgal 3 , B = 24 


( 67 ) 


- 26 - 


ERROR (cm) 


Figure 8. Total RMS Errors for 1° x 1° Mean Undulation Differences, 


Assuming Perfect Gravity Anomaly Data Within a Cap of 
Radius and Potential Coefficients to Degree 180 Estimated 
from 1° x 1° Gravity Anomalies Accurate to . 
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The parameters for the model (66), as given by Jekeli (1978), are 

a x = 18.3906 mgal 2 , cr 1 = .9943667 A = 100 

(X 3 =658.6132 mgal 2 , a a = .9048949 B = 20 

This model refers to the sphere of radius R (cf. p. 9 ). Comparing the two 
models, in general, it is clear that for large n , the first term of model (66) 
provides the dominant contribution; and therefore, the decrease of c n in (66) 

is less rapid (for 100 < n < 600) than in (67). The implication is that the 

truncation error, which is essentially the sum of the high-degree error vari- 
ances, is greater- when model (66) is used instead of model (67). This conclusion 
is verified by the values of Table 3 which are computed via equation (37) ( Method 
A ) with n = 180 and $ 0 = 0° . Similar results would evidently be obtained if 
other cap radii and equation (36) (Method B) were used. 

Table 3. Truncation Errors of Undulation Differences (00 = 0°, n = 180, 
Method A ) based on two different c n -Models. 


0 pc; (km) 

. — — . — . — - 

m 2 ' ( cm) 
model (67) 

m^ ( cm) 
model (66) 

200 

63.6 

88.3 

600 

66.6 

93.2 

1000 

67.1 

94.0 

1400 

67.1 

93.9 

1800 

69.9 

93.6 


It cannot be stated categorically that one of the above models is superior to 
the other since the characteristics of the gravity field on which they differ are 
not yet established with sufficient confidence. Therefore, the considerable 
differences between the last two columns in Table 3 suggest still some uncer- 
tainty in the estimation of the truncation errors. 

10. Other Methods 

By no means is the method based on Molodenskii's truncation theory the only 
approach to undulations from gravity disturbances. Collocation ( Moritz 1972) is 
a method as amenable in this case as when gravity anomalies constitute the gravity 
data. Although this method is not discussed further (Christodoulidis (1976) shows 
that the corresponding error analysis is quite difficult for large amounts of data), 
the cross -covariance function between gravity disturbances and geoid undulations 
is given here for the sake of completeness (cf. equation (25)): 

oov (Np, 6a) = JEL V iti. c, (-Sl)" +1 p n ( cos (lip, ) (68) 

>r Q n A,(n -!) \r P r Q / 
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The sum starts with n = 2 under the assumption that the zero- and first- 
degree potential harmonics are absent; also c n refers to a sphere of radius R. 

An error analysis was conducted by Sjoberg (1979), in which the geoid 
undulation is assumed to be computed from global mean gravity data. The lack 
of more detailed data than mean anomalies constitutes a type of truncation 
error; while erroneous data implies a commission error. The global mean of 
these errors in this case can be written in terms of the smoothing factors 8 n . 
The commission error also relies on an error model similar to equation (53). 
Both error estimates are infinite series and are also developed for undulation 
differences by Sjoberg (1979). On the basis of a 1° x 1° mean anomaly field 
with an accuracy of 1 mgal, computations using these formulas give a total 
error in AN of 55 cm ( 0 pq = 2000 km). Assuming that the lack of detail beyond 
1° x 1° mean anomalies is equivalent to a truncation error with n = 180, the 
curve for tyo- 0, e^ g = 1 mgal in Figure 6 shows an error in AN of 68 cm 
(0 pq = 2000 km). The difference in these results (~20%) is not excessive in 
view of the totally dissimilar methods of analysis. 

It was noted earlier ( section 3 ) that the Molodenskii truncation theory 
provides for the series of a part of Stokes' integral to converge faster than the 
series expansion of the entire integral. Molodenskii, et al. (1962) discuss a 
further improvement in this convergence. A very elegant type of generalization 
of their principle was investigated by Colombo (1977). The concept of this meth- 
od rests on the simple fact that a change in the Stokes’ function ( the kernel of 
the integral) naturally induces a change in the errors as embodied in part by 
the (now also changed) Q n coefficients. Colombo parameterized the kernel 
and achieved, to some extent, the "best" Stokes' function by minimizing the 
errors in a least squares sense. These promising techniques can obviously be 
applied to the integral of gravity disturbances without any difficulty. 


11. Summary and Conclusions 

It was the intent of this report to subject the geoid undulation differences 
as they are determined by a surface integral of gravity disturbances to an error 
analysis in the frequency domain. The principles of the Molodenskii truncation 
theory are applied to the extent that the undulation N is decomposed into an 
integral over a spherical cap centered at the computation point and an infinite 
series of spherical harmonics which essentially represent the contribution to 
N from the gravity field outside this cap. A study of the data requirements as 
they pertain to the acquisition of gravity disturbances within the cap was omitted 
under the assumption that they do not differ significantly from the correspon- 
ding requirements of gravity anomaly data, as examined by Christodoulidis (1976). 

The errors of primary interest, therefore, were those due to the neglect 
of the high-degree potential harmonic coefficients, as well as the errors in the 
available coefficients. A new model for the standard deviations of the coefficients 
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was introduced to better simulate their accuracy when they are determined 
from terrestrial gravity data. The final analysis was based on an assumed set 
of potential coefficients to degree 180 and on various cap radii between 0° and 
30 & . Furthermore, the analysis was expanded to include investigations with 
respect to the conventional formulation of point and mean undulations; that is, 
in terms of gravity anomalies within the cap plus potential coefficients. 

Figures 2 through 5 demonstrate that as the distance between the points of 
computation increases, so do the commission errors ( due to erroneous potential 
coefficients) for all values of the cap radius, as well as h, the maximum degree 
of the coefficients. However, for large values of n ( e.g. n = 180 ), the trun- 
cation error is essentially independent of the separation between the points. 
Naturally, both the commission and truncation errors increase with a decrease 
in the cap size because the undulations must then rely more heavily on the po- 
tential coefficients. 

A comparison of errors derived through Methods A and B ( see (56)) showed 
that Method B, for the most part, is not inferior to Method A, and in some cases 
it may yield as much as a 30% improvement in accuracy ( cf. the curves in Fig- 
ure 7: for 0 P q = 1500 km, = 10 mgal (1 mgal), Method A gives 48 cm ( 8 cm) 
and Method B gives 33 cm ( 5 cm) ). For small spherical caps ( 0 o < 2°), the two 
methods give practically equivalent error estimates ( Figure 7: the error for 
0 O - 2°, 0p Q = 1500 km, e^ s = 10 mgal ( lmgal) is about 75 cm( 13 cm); with 
a maximum error (at 0 O - 0°) in undulation differences of 151 cm ( 68 cm )). 

The confidence that can be placed on either method taken individually depends on 
how realistic the models are for the degree variances of gravity anomalies and 
potential coefficient errors. For example, the truncation error can vary 30% to 
40% by changing the degree variance model c n ( see Table 3). ( However this 
may be of little consequence to the total error when n is large, for then the 
truncation error is often smaller than the commission error.) 

As mentioned in section 1, the possibility of determining oceanic geoid 
undulations from an integral of gravity disturbances, thus dispensing with reduc- 
tions to the geoid, prompted the investigation of Method B. In contrast. Method 
A is predicated on the availability of gravity anomalies reduced to the geoid, this 
being a feasible hypothesis if the sea surface topography is known beforehand, 
perhaps through an oceanographic model ( or if the error caused by neglecting it 
can be tolerated). For small caps, when both methods of computing undulations 
render similar errors, the route to follow may be decided by weighing the prob- 
lem of acquiring such sea surface topography against the need for altimeter 
measurements which are required for the computation of gravity disturbances. 

The latter requirement demands considerations with respect to the accuracy of 
the radial positional component of a satellite -borne altimeter ( and the accuracy 
of the altimetry itself). For example, in view of the radial component of the 
gradient of normal gravity (% . 31 mgal/m), the sea surface height above the 
ellipsoid should be known ( from altimetry, see Figure 1) to better than 3 m in 
order to derive gravity disturbances with an accuracy greater than 1 mgal ( this 
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being required for data within the cap if the potential coefficients for the outer 
zone contribution are determined from mean anomalies accurate to 1 mgal; 
see section 7). 

For larger caps ( = 5° to 30°), the apparent improvement in the trunca- 

tion and commission errors according to Method B may constitute an additional 
factor to consider when choosing the particular method of computing undulations. 
Quite significant reductions in the error, as compared to Method A, are found 
for large caps when the zero- and first-degree harmonics of the gravity distur- 
bances do not exist ( or are known) (and this may be the more appropriate prem- 
ise on which to base comparisons to Method A). Appendix C includes an abbre- 
viated analysis of errors derived through the modified kernel function S*($) 

( see section 3 ). It is shown there that for ^ 0 = 30° and small separations $>q , 
the reduction in error ( from Method A to Method B, modified) is 75% to 90%, 
and for $ 0 = 10°., it is as much as 50% . For larger separations ( i/) P q> 1000 km), 
the difference between Method A and Method B (modified) is substantial (about 
50%) only if the truncation error is the predominant error source (e.g. if = 

1 mgal). On the other hand, as increases, the difference between Method 
B and Method B (modified) becomes less conspicuous. In the comparisons above, 
it is assumed that the gravity data within the cap is perfect and continuous. The 
reader is referred to the table of Appendix C for specific values of the estimated 
RMS errors. - 

Finally, in regard to mean undulation differences, it was seen that the 
estimated errors are substantially reduced only when the truncation error is 
the dominating constituent. In the latter case, a typical reduction of 50% ( from 
68 cm to 31 cm, for = 1 mgal) is attained when potential coefficients to 
degree 180 are used to compute l°x 1° mean undulation differences instead of 
corresponding point values. 

Taking all numerical comparisons into consideration, it is apparent that 
for small caps (0 O < 10°) the stringency of the data requirements as estimated 
by Chris todoulidis (1976) for a 10 cm relative oceanic geoid has not diminished 
significantly in the case when undulations are determined on the basis of gravity 
disturbances. For larger caps, some improvement in the truncation and com- 
mission errors is indicated, but whether sufficiently dense and accurate gravity 
data can be amassed within large caps ( say, with radius 0 O = 20° or 30° ) is 
questionable at present. 
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Appendix A 


The integration of 

rZ. 

o) = ,S(0)P n (cos 0) sin 


is achieved by utilizing the recurrence relations of Legendre polynomials. Let 
P n (x) (or simply P n ) denote the n to -degree Legendre polynomial and let P n '(x) 
(or simply P n ') denote its derivative with respect to the argument x. The fol- 
lowing recurrence formulas are well known: 


0 d0 , 0o * 0 


Qn(0 


(n + 1) P n+1 

+ nP n _i 

= ( 2 n + 1 ) x P n 

(A.l) 

(l-x S )P n r 

= n (P n _! 

- xP n ) 

(A. 2) 

( 2n + l ) P n 

= PnVl * 

Pn— l 

(A. 3) 


To facilitate the integration, let t = cos0 and x = cos 0 O . Then the function 
S(f) ( see equation (4)) becomes 


and 


S(t) = 
Qn(x) = 


1-t 

i»X / ip 


&n 


(WA)’ tf61 


(A.4) 


rt X 


1 -t 


P n (t)dt - J^p n (t)072(l+y^jdt, x^ 1 (A. 5) 


In all further derivations, it will be understood that x ^ 1. For n> 0, 

... 1 T X Pn(t)dt 1 P X P„(t) 

let L „<x) = -$=±=r , K n (x> = iTf JjTTLp dt 

andfornil, I„(x) = [ P.(t) dt = — i— ( P n+1 (x) - P,.! ( x) ) 

J i 2 n + 1 

The last equality follows from equation (A. 3) and the fact that 
P n ( -1) = (-l) n , for all n :> 0 
Hagiwara (1976) shows that 

K n + 1 - 2K n + K n _! = I n , n > 1 

1 / / O \ 

with K 0 


(A. 6) 


(A. 7) 


and* 


2 V 1 

I, 


/ 2 ( 1 - x) 

2 \ / 

- K 0 - (^1 - 

1 


1 -x 


Ki 


1 — x 


JL 

jn /2 ( 1 - xj ~ t K n + l~ K n-l] 2 n + 1 V n ~ 1 

* This differs from the misprinted version in Hagiwara's paper. 


1 (A. 8) 
(A. 9) 
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The first integral in (A. 5) is 
, x rr 


J 2 

P n( fc ) dt = ^ L n ( X) 

Let the second integral in (A. 5) be 

^ 

: n (x) = - | p l( t) & n (l •• 1 2 


Rn 


1-t 


dt 


(A. 10) 


(A. 11) 


Substituting the recurrence formula (A.l) gives 

r„<x) =— R 0 - a( x ) -^ij i tp n . l(t) ^( 1 + ^ir) dt (A#12) 

Now, putting (A. 2) in the integral of (A. 12) results in 


n-1 


2 n - 1 f’ x j-. 


R„ (X ) - -^h^x, -“j (iVrr) + 


n 


2 \ 


dt 


= R "- s(x) + '^f^iT J_i (1 ' t2)P " Li(t)fe ( 1+ '/r?7) dt < A - 13 > 

The integral in (A. 13) is integrated by parts. To this end, let 

* = d-t 2 ,c».( 1+ y^i) 

After some simplifications, it is found that 

dw = ["/ I r +1 - 2t0 ”( 1+ VS)] dt 

Also, if dv = P n Lx (t) dt, then v = P n _ x ( t) and the integral in (A. 13) becomes 


f 

J -1 


w dv = w v 


|X r»x 

-1 J-l 


v dw 


= <l-x 2 )^(l + 7^)p n . l(X ) + iJ/ i -TTp n . l(t)dt + 

- p nrl( 4 ) dt +2 P n .,(x) + J-S-'j dt 

Substituting this back into (A. 13) yields 


R n (x) - R n -a( x ) + — ~ ( 1 -x 2 ) 8m (l + 


n(n-l) 


ra p - i(x ) + (a. i4) 


+ 


2n ~ 1 A 2n-l 2n-l 

t/2(n-l)n A -i( x ) “ n ( n - 1 ) Tn ~ 1 ^ x ^ + 2 ^n -l) B n-i( x ) 


T n-i (X) is given by (A. 6), and 

A «-l(X) = /_* /Pfp tl (t) dt. 
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B n _ 1 (x) = J X tp^^t) Qn(l+ JjL ^ dt 

A n _ i ( x) is computed by inserting (A. 3) and integrating by parts: 

A„-i(x) = [/T^T(P n (x) - P n _ a (x))-/T(I>(-l)-P n _ 3 (-l)) + 

+ ( L n ( x) - L n _s( x))] (A. 15) 

The second term in (A. 15) is zero in view of (A. 7) and L n (x) is given by (A. 6). 
Finally, B n _ a ( x) is obtained by substituting (A. 1): 

B n— 1 ( x ) = - ~~ R„ ( x ) - “'Tf R„- 2 (x) (A. 16) 

Also, from (A. 6) we find 

In-l(X) = -^Tl(P n ( X) - Pn-2<X) ) (A . 17) 

Putting (A. 15), (A. 16), and (A. 17) into (A. 14), it is readily found that 

= ^r„- 2 <x )+ ^i_ ij(Mx) + 

+ n(n-l)/r < /Tr ^~- /2~)(P»(x)- P,- S (x)) + 

+ ^fj^( 1+ ,/i?i)( 1 - x2 ) p ->( x > < A ‘ 18 > 

The remaining problem is to find a suitable recursive relation for the L n ( x) 
functions. This can be done with some juggling of equations (A. 8) and (A. 9). 

For convenience, let u = 1 / /2 ( 1 -x) , then 

(A. 8) becomes K n + 1 - 2K n + K „_ l = -u I n (A. 19) 

and (A . 9) becomes L n = u I n - ( K„ +1 - K n _ x ) (A . 20) 

2n+ 1 

Solving (A . 19) for K n+1 and substituting into (A. 20) yields 
(2n + l) L n ( 2n + 2 ) u I n - 2K n + 2K n _ x 

Also, ( 2n + 3 ) L n+1 - (2n + 4) u I n+1 - 2K n+1 + 2K n 
Adding these last two identities and substituting (A. 20) gives 

(2n + 3) L„ +1 = 2u [ (n+2 ) I n+1 - n IJ + (2n + l) L n 

Or, finally 
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(2n+l)L„(x) = 7 [(n + 1) I n (x) -(n-1) I n _ 1 (x)] + 

+ ( 2 n - 1 ) L n _!(x), n ^ 1 

with L 0 (x) = 2 - 2 (A. 21 ) 

Lq(x) is easily determined from (A. 6). The starting functions R 0 (x) and R x (x) 
are found by integrating (A. 11) directly, changing the variable t to 

t = 1 - 2 z 3 , z ^ 0 

and integrating by parts. Only the final results are stated below. 

The recursive formulation of Q n may be summarized as follows: 

For x = cos 0 O , (6 0 ^ 0 : 

Q n (x ) = 2 L n ( x ) + R n (x) , n £ 0 

where (2n + l)L n (x) = [(n + l)I„(x) - (n- 1 ) I^x) J + 

+ ( 2 n - 1 ) L n _!(x), n ^ 1 

with L 0 = 2-2 

and : " (x) “ nhi i p »+i( x > - p »-i( x >]* n = 1 



and (n + 1) R n ( x ) = i(n-l)(n- 2 )R„_ a (x) + L„(x) -L„- a (x)] + 

2 n-lr /l -x ,1 

+ “"ITl V ~ " 1 J In - 1 ( x ) + 

+ ^ S ^< 1-x2)Pn -^ x > 0?I ( 1+ Vi?7)* n s 2 

with R 0 (x) = + + + 

and R l( x) = ( 1 _i|i)[(i-x,e, ! (l + ^ll)-l] + |(l- 


— X 
2 


1 — x X — x 


and (n + 1) P B)1 (x) +nP H (x) = ( 2n + l ) x P n ( x), n ^ 2 

with P 0 (x) = 1 
and P x ( x ) = x 
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Appendix B. Concerning the factor cr R in equation (46) 


The extension of the disturbing potential into the space outside a sphere 
of radius R B is given by (cf. equation (17)) 

CO t! 

T (r, 9 , A) " —■ V £ ( C na cos mA + S^sinmA ) P w (cos 8) 

nsd a rsQ 


, V T 

} in Vr / 


1H-1 


/_! 

n 0 


(B.l) 


where C„«> ST nll refer to this sphere ( that; is, R is replaced by Bo inequation 
(17)). On the sphere of radius R > R n , it is 

T(R,e,X> -f T. (f^ « f T. (B.2) 


II “0 ' ft ™ 0 

If N 3 is computed on the sphere of radius R, then according to equation (16) it 

IcM 

if 


is 


No 


R 
2 y 


Z i 

* 


y 


(B. 3) 


n~0 


where 5g u on this sphere is obtained from (18) and (B.2) : 


n + 1 
R T » 


n-M 

R 


T, 


Rb\ 
•« \R / 


n +1 


(B.4) 


Consequently, the substitution of (B.4) and T n of (B.l) into (B.3) leads to 

«■» „ n 

V*' P V* cn 

N a « ) (n + 1) Q„ r~) } (C nn cosmA +S 0 .sinmA) P M (cas 8) 

M n ~0 K ^ 

The corresponding formulation of the error differs from equation (38) by the 
factors ( Rr,/ R f ; this implies that the degree variances £ n as given by (40) 
should be modified by the factor cr n -• (Rtf/ B a ) u . 


Appendix C . Truncation Coefficients Implied by the Kernel Function 
S*(0) « S(0)-l-f cos ip 


If the series expansion of the gravity disturbance function 6g- ( equation (14)) 
contains no zero- and first-degree terms, then the integral equation (5) becomes 

N “ 4^7 JJ l *•& ^ p * <eo8 *> *•*> 

Or n = a f -3 

who ix' the series expansion for S((!)) (equation (4)) has been used and l^(cos (l)), 
& ■- 0, 1 are excluded because they are orthogonal to the harmonic funct ions of 
Sg. Equivalently, equation (C. 1) is 


N " 4^ J J 5 K C 6 ' d0 


( 0 . 2 ) 
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where S *(ip) = S ( 0 ) — 1 — §cos Ip. This modified kernel function defines the 
following truncation coefficients (cf. equation (11)) 


Qn* ~ f S*(0) P„(cos^) sin0 d0 
J 0o 


= Qn" 


TT 


00 


( 1 + §• cos 0) P n (cos0) sin 0 dip 


(C.3) 

(C.4) 


It is noted that for 0 O = °* Q n *=Q n , ns 2; also Q 0 * and Q x * are not equal 
to zero unless ip Q = 0, but this is immaterial in the error formulas since 5g 0 = 

0 = 6gi. 


The integral in equation (C.4) is evaluated by utilizing the results of 
Appendix A ( let t = cos0, x = cos lp Q ). 


.rr 

Pn(cos Ip) sill ip dip = I n (x) 

00 

, TT «X 

cos 0P n (cos ip) sin 0 dip = tp n (t)dt 
0o -1 


(C.5) 

(C.6) 


The integral in (C.6) is directly integrable when the recursion formula (A.l) is 
introduced: 

X X 

tp n (t)dt= — - f ((n + 1) P n + 1 (t) +nP n _ 1 (t))dt 
J_X 2n + l J-x 


2n + l + 2n + l (Q*^) 

Substituting (C. 5) and (C. 7) into (C.4), 

Qn* = Qn " [ In + *T ( (B+l) Ul + nl n _ x ) ] 

4n + 2 

where Q n and I n are given on page 38. 

The equation for the global error estimates m x and m 2 of the outer zone 
undulation differences are obtained by replacing Q n with Q n * in equations (46) 
and (36). The following table indicates the type of corresponding reductions that 
can be expected in the error estimates. 


- 40 - 



Total ( RMS ) errors (cm) in undulation differences assuming perfect and continuous 
gravity data within a cap of radius ip 0 ( n = 180, Method B, parenthetical values 
are obtained via Method A ). 




= 1 mgal 

<r 4s = 10 mgal 


tppq [km] 

using Q n (Q n ) 

using Q n * 

using Q n (Qn) 

using Q n * 


100 

4.8 (7.4) 

3.5 

8.2 (12.4) 

mm 

-e- 

o 

li 

M 

O 

0 

mm 

5.0 (7.6) 

3.8 

27.3 (39.5) 




5.7 (8.6) 

4.7 

37.9 (54.8) 

m 


mm 

2.1 (1.7) 

0.1 

3.6 ( 3.1) 

0.8 

tf>0 = 30° 

BSH 

2.3 (2.1) 


12.9 (14.5) 

8.3 


EM 

2.6 (2.9) 

1.6 

18.7 (25.0) 

16.2 


Generally, for most n, Q n * is less in magnitude than Q n , and the absolute 
differences between these coefficients increase substantially as the radius iji 0 
of the cap changes from 10° to 30°. The truncation error (n = 180) for tf) 0 = 30° 
is practically nonexistent (« 0.02 cm, compared to 1.9 cm when Q n is used); 
while for 0 O = 10°, the truncation error of 3 cm to 4 cm is closer to the 4 cm to 
5 cm that is obtained with Q n . This accounts for the considerable decreases in 
error for the larger cap, particularly when 0 PQ is small as the commission er- 
ror is then minimal. The tabulated values which are enclosed in parentheses 
correspond to error estimates based on Method A (using the coefficients Q n ). 
Note that both Q n and Q n * are derived from kernel functions (S ( 0) and S*( 0), 
respectively) whose zero- and first-degree harmonic coefficients are zero. 

We see then the noticeable improvement (up to 50% for 0 O = 10°) in the 
accuracy of geoid undulation differences when determined on the basis of perfect 
and continuous gravity disturbance data within a cap of radius 0 O - 30° , and 
when the modified kernel function S *(ip) is employed. This naturally leads to 
further possible improvements with more discriminating selections of the kernel 
function, as alluded to in section 10. 













